function derivative = regressDeriv(x, y, window, weights)
  % Function derivative estimation via simple linear regression
  % Given (in)dependent variables (x) resp. y, estimate the derivative y' of the function
  % y(x) at x_i, y'_i, as the slope of the least-square-error linear fit of the points
  % (x_-window/2, y_-window/2) ... (x_window/2, y_window/2).  weights (optional) is a
  % numel(x)-numel vector of weightings
  %
  % Author: David Goldsmith, Wash. State Dept. of Ecology, dgol461@ecy.wa.gov
  % Release date: 9/15/2011
  
  % Get numel and check data consistency
    mnx = size(x);
    mny = size(y);
      % x & y must be 1D & have more than one element
    if (numel(mnx)>2) || (numel(mny)>2) || ...
       (sum(mnx==1)~=1) || (sum(mny==1)~=1)       
        derivative = ['Sorry, one or both of the inputs to regressDeriv '...
                      'is the wrong size!'];
        return
    end
      % x & y must be "column vectors"
    if mnx(2)~=1, x = x'; end
    if mny(2)~=1, y = y'; end

    npts = numel(window);
    
  % If necessary, set weights
    if nargin < 4
        weights = ones(size(x)); % weights has to have same orientation as x...
        weights = weights(1:npts); % ...but same number of points as window 
    elseif isempty(weights)
        if rem(npts,2) % if npts is odd
            weights = [1:ceil(npts/2) floor(npts/2):-1:1];
        else % npts is even
            weights = [1:npts/2 npts/2:-1:1];
        end
    end
    
  % Now regress down the line
    derivative = nans(size(x));
    for i=1:size(x,1)
        try
            range = window + i;
            X = [(weights .* x(range)) weights];
            Y = [weights .* y(range)];
            regress = X \ Y;
            derivative(i) = regress(1);
        catch
            continue
        end
    end
end